Assignment   a3_density plots				    
By           Daniel B. Carr 
Copyright    2005,2006,2007 permission granted for class use
                                
Goals        Understand kernel-based density estimates
             Use data to estimate densities 
             Investigate parameters of theoretical density function 

Assumptions  Observations are random sample from
             a continuous distribution. 

Exercises    1 Constructing a kernel density estimate
               example with 5 points
             2 R density plots for data
             3 Theoretical Gamma density plots 
             4 Experimenting with Beta shape parameters

Read         Comments related to Density Estimate will soon be available
             on the class web site.  
		
Due          1. Fix part 1.11 and submit plot
             2. Plot
             3. Plot
                Look for directions near the end of
                the script. This is set up for you
                to use the locator() function to put
                labels at mouse click locations in the plot.                      
             4. Plot and answers to four questions.

1. Kernel Density Estimation

Here the goal is to provide an intuitive understanding of
kernel density estimates.  

There are many ways to obtain densities and sophisticate
ways to obtain smoothing parameters for kernel density estimates.   
For those interested in learning more I recommend
David Scott's book on Multivariate Density Estimation,
John Wiley & Sons as an introduction to a vast statistical
literature.  
    
The ideas behind kernel density estimation are simple.
  
1) We treat each observation as surrogate for unobserved data
   in a local neighborhood of the observation 
2) We attach a local density function to each observation.
   This local density function indicates our guess at the
   relative likelihood for the local unobserved data  
3) We construct a composite density function as the average
   of the all the local density functions.  At any location
   on the real line we average the specific values from the 
   local density functions.  

1.1 Graphics colors, device and parameter___________________

## Run BE SURE to FIX the ??'s and UNCOMMENT the text command
##     just before the ##End  

windows(width=5,height=6)
gridY = c(.0,.05,.1)
gridLab = c(0,.05,.1)
cex = .6

#1.2 Panel layout _______________________________________

n = 5   # number of samples and layout parameter
layout1 = panelLayout(nrow=n+1,ncol=1,
          topMar=.5,bottomMar=.3,
	    leftMar=0,rightMar=.3,
          rowSep=c(rep(0,n),.20,0))

#1.3 Observations locations_________________________

# Sorted sample values were obtained using the line below.  
# obsLoc = sort(rnorm(n,mean=100,sd=10))

# Some samples are undesirable for illustration purposes so the
# command repeats until the value below were obtained 

obsLoc = c(77.24799, 97.59442,106.23001,106.63648,110.65806)

# Rather than explicit retype the values as above, an alternative
# is to same the random number seed just before getting the sample.
# Reinstating the seed would case the same sample to be generated. 

#1.4 Obtain local relative likelihoods values

# We will use a normal kernel for the local likelihoods in our density estimation

cumProb = ppoints(100)			         # cumulative probabilities
sd = 3 				               # The standard deviation of the normal kernel
kernelQuants = qnorm(cumProb,mean=0, sd=sd) # quantiles
kernelDensity = dnorm(kernelQuants,sd=sd)    # densities or ordinates

#1.5 X-axis scaling____________________________________________________

obsRx = range(obsLoc)
obsRxExtend = obsRx + range(kernelQuants)	  # Extend: add room for extreme kernels
# pad the scale for plotting
rx = mean(obsRxExtend) + 1.04*diff(obsRxExtend)*c(-.5,.5) 

#1.6 Locations for drawing vertical stripes__________________________

nstripe = 50 # The number of location estimates
stripe.rx = obsRxExtend
#stripe.rx = rx + .01*diff(rx)*c(-1,1)	  # the range for vertical stripes
stripe.loc = seq(stripe.rx[1],stripe.rx[2],length=nstripe)

#1.7 Y-axis scaling_________________________________________________________

ry = range(kernelDensity)    # range of heights for the kernel 
kmin = ry[1]		     # minimum kernel height represented
ry[1] = 0		           # set min to 0 and subtract from kernel
ry[2] = 1.04*diff(ry)	     # expand the range by 4%
ry[2] = .14		           # alter for nice tick labels
 
#1.8 Draw kernels and stripes________________________________________________
i = 1
gridLab = c(0,.05,.1)
for (i in 5:1){
   panelSelect(layout1,i,1)
   panelScale(rx,ry)
   panelFill()
   panelGrid(y=gridLab[-1],col="white",lwd=2)

   # draw black bell curve polygon
   polygon(obsLoc[i]+kernelQuants,kernelDensity-kmin,density=-1,
   col='blue') 

   # stripes y for all locations
   stripe.y = dnorm(stripe.loc,mean=obsLoc[i],sd=sd)

   # values below kmin are outside the local bell curve 
   good = stripe.y > kmin
   stripe.x = stripe.loc[good]
   stripe.y = stripe.y[good]-kmin  # a small adjustment 
   stripe.n = length(stripe.x)
   polygon(obsLoc[i]+kernelQuants,kernelDensity-kmin,density=-1,
      col="black")
   segments(stripe.x,rep(0,stripe.n),stripe.x,stripe.y,col="white")
   polygon(obsLoc[i]+kernelQuants,kernelDensity-kmin,density=0)

   points(obsLoc[i],0+.055*diff(ry),pch=17,col="white",cex=1.5)
   mgp = c(2.5,.4,0)
   axis(side=4,at=gridLab,labels=substring(format(gridLab),2,4),
       mgp=mgp,tck=0,cex=cex,las=2)
   panelOutline()
}

#1.9 Averaging ordinates (stripe heights)_______________

stripe.ave = rep(0,length(stripe.loc))
for(i in 1:n) stripe.ave = stripe.ave+
	           dnorm(stripe.loc,mean=obsLoc[i],sd=sd)
stripe.ave = stripe.ave/n

#1.10 Draw the density estimates_____________________________

panelSelect(layout1,n+1,1)
panelScale(rx,ry)
panelFill()
panelGrid(y=gridLab[-1],col="white",lwd=2)
polygon(c(stripe.loc[1],stripe.loc,stripe.loc[nstripe]),
        c(0,stripe.ave,0),density=-1,col="black")

segments(stripe.loc,rep(0,nstripe),stripe.loc,stripe.ave,col="white")
polygon(c(stripe.loc[1],stripe.loc,stripe.loc[nstripe]),
        c(0,stripe.ave,0),density=0)

panelOutline()

axis(side=4,at=gridLab,mgp=mgp,labels=substring(format(gridLab),2,4),
     tck=0,cex=cex,las=2)
axis(side=1,at=c(70,80,90,100,110),tck=-.04,mgp=c(2,.3,0),cex=cex)

#1.11 Add a title in the top margin panel________________

panelSelect(layout1,mar='top')
panelScale(c(0,1),c(0,1))

# !!!Uncomment and fill arguments designated with question marks.
#    Center the text center
# text(?,?,'Density Estimate Construction',adj=?,cex=?)

##End 

2. R density() plots for data_______________________________________

R provides a choice of kernels for kernel density
density estimation and provides some algorithms
for selecting a smoothing parameter.  There
is a lot written on these topics.   

Conventional wisdom that the precise choice of kernel
does not matter that match.  Some theory suggests
that the Epanechnikov kernel has special merit.  
The default "gaussian" is fine for our purposes.

Conventional wisdom is that the important choice
is the smoothing parameter.  The smoothing parameter
for kernel estimates can be thought
of the standard deviation of kernels that are averaged 
in producing estimated density function as in example 1.
The smooth parameter often called the bandwidth
R provides several methods prefixed with "bw." (bandwidth)
for picking the bandwidth for a gaussian kernel based 
on the data, x.  Here we use methods bw.SJ()  

##Run
?bw.SJ      # See the documentation in R  
##End 
 
2.1 Digression us dotchart() look at data used in the kernel

The exercise uses precipitation data the comes with R.
##Run
?precip
##End  

##Run
windows(width=8,height=10)
require(graphics)
dotchart(sort(precip),
    main="Average Annual Precipitation for Cities, Period ??",
    cex=.7)
title(sub="Precipitation in Inches")
##End

2.2  A density plot  

##Run

# Find a decent smoothing parameter for the datat
bw = bw.SJ(precip) ## sensible automatic choice


# Try several different kernels.  
# The chosen smoothing parameter may be poor for some of them
kernels = eval(formals(density.default)$kernel)

windows(width=8,height=10)
panels = panelLayout(nrow=7,ncol=1,
          leftMar=1.5, rightMar=1,
          topMar=.4, bottomMar=1,
          rowSep=.1)

# Default gaussian kernel, used to find a range
ans = density(precip,bw=bw) 
rx = range(ans$x)
ry = range(ans$y)
ry[1] = 0
ry[2] = ry[2]*1.10   # pad by 10%

for (i in 1:7){
   panelSelect(panels,i,1)
   panelScale(rx,ry)
   panelFill()
   axis(side=1,tck=1,labels=F,col='white',las=2,lwd=2)
   axis(side=4,tck=0,mgp=c(2,.3,0),las=2)
   axis(side=2,tck=1,labels=F,col="white",las=2,lwd=2)

   lines(density(precip, bw = bw, kern =  kernels[i]),lwd=2)

   mtext(side=2,line=1,las=2,kernels[i],adj=1)
   panelOutline()
   if(i==1)mtext(side=3,line=1.5,"Densities Based on Seven Different Kernels",cex=1.2,adj=.5)
   if(i==7){
      axis(side=1,tck=0,mgp=c(2,.3,0))
      mtext(side=1,line=1.5,"Average Annual Precipitation in Inches")
   }
}  

##End


3. Overlay of Gamma density plots___________________________________________

##Run
cumProb = ppoints(50)
q1 = qgamma(cumProb,shape=1)
d1 = dgamma(q1,shape=1)

q5 = qgamma(cumProb,sh=5)
d5 = dgamma(q5,sh=5)

q15 = qgamma(cumProb,sh=15)
d15 = dgamma(q15,sh=15)

q35 = qgamma(cumProb,sh=35)
d35 = dgamma(q35,sh=35)

rx = range(q1,q5,q15,q35)
ry = range(d1,d5,d15,d35)

windows()
plot(rx,ry,type='n',axes=F,
     main='Theoretical Density Plots',
     ylab='Density',
     xlab='Gamma: Shapes=1, 5, 15, 35')
panelFill()
axis(side=1,tck=1,labels=F,col="white",lwd=2)
axis(side=2,tck=1,labels=F,col="white",lwd=2)
axis(side=1,tck=0,mgp=c(2,.3,0))
axis(side=2,tck=0,mgp=c(2,.4,0),las=2)
lines(q1,d1,lwd=2,col="blue")
lines(q5,d5,lwd=4,col="cyan")
lines(q15,d15,lwd=5,col="magenta")
lines(q35,d35,lwd=6,col="yellow")
panelOutline()

#! move the plot were it is convenient to 
#  mouse click on locations in the plot 
#  

locs = locator(4)  # Click on four point near the
                   # respective peaks left to right   

text(locs$x,locs$y,paste('Shape = ',c(1,5,15,35),
     sep=''),adj=0)

##End

3.1 Comment:  Gamma Distribution and
          the Central Limit Theorem for Sums_________________

The gamma distribution characterizes the total waiting
time for n (shape parameter) independent events from a fix
rate Poisson distribution.  

For a single event (shape parameter 1) the waiting time
follows an exponential distribution

The waiting time for n events is the sum of waiting
times for n first time events.  

The commonly stated Central Limit thereorm applies
to the average of independent identically distributed
random variables.  The sum of variables is the average
multiplied by the scale factor n,  the number of items
averaged.  The sum will also have a normal distribution
as the sample size n gets large.  

Observe the bell like shape of when the Gamma distribution
has shape parameter 35.

4. The Beta Distribution_____________________________________________

The support for the standard beta distribution is the inteval
[0 1].  This can be adapted to cover different intervals
by adding two parameters to set the minimum and maximum values. 

The Beta distribution on [0 1] has two shape parameters.
These let the beta distribution take many forms. 
The key part of the density is
   x^(shape1-1) * (1-x)^(shape2-1).

The special case, shape1 = shape2 = 1, has 
exponents of zero.  The density is 1 over the
interval [0 1] so it is  uniform distribution. 

As a reminder, the four
basic distribution functions are  

density:                dbeta(p, shape1, shape2) 
cumulative probability: pbeta(q, shape1, shape2) 
quantiles:              qbeta(p, shape1, shape2) 
random numbers:         rbeta(n, shape1, shape2) 

The idea here is to do a little experiment
to see how the shape parameters change
the density.

##Run
windows(width=7.5,height=7)
shape =  c(.5,1,2)   # Our experimental values

# set the y scale for comparabilty
cumProbs = ppoints(50)
xGrid = c(.2,.4,.6,.8)
ry = c(1,1)  
for (i in 1:3){
for (j in 1:3){
  betaDensity = dbeta(cumProbs,shape1=shape[i],shape2=shape[j])
  ry = range(ry,betaDensity)
}}
ry[2] = ry[2]*1.05
shapeText = format(shape)

pan = panelLayout(nrow=3,ncol=3,
      leftMar=.5,rowSep=0,colSep=0)

for(i in 1:3){
for(j in 1:3){
   panelSelect(pan,i,j)
   panelScale(c(0,1),ry)
   panelFill()
   axis(side=1,tck=1,labels=F,col="white")
   y=dbeta(cumProbs,shape1=shape[i],shape2=shape[j])
   polygon(c(cumProbs,1,0),c(y,0,0),density=-1,border=F,col="skyblue")
   lines(cumProbs,y,col="Blue",lwd=3)
   panelOutline(lwd=2)
   if(j==1)mtext(side=2,line=.6,shapeText[i],las=2)
   if(i==1)mtext(side=3,line=.4,shapeText[j])
   if(i==3)axis(side=1,at=xGrid,mgp = c(2,.1,0),tck=-.04) 
}}

panFull = panelLayout(nrow=1,ncol=1,
                      leftMar=.5)
panelSelect(panFull,1,1)
panelScale()
mtext(side=3,line = 3.3,"Beta Densities",cex=1.2)
mtext(side=3,line=1.6,'Shape2 Parameters')
mtext(side=2,line= 2.1,'Shape1 Parameters',lwd=2)
##End

Answer the four questions based on the plot  

From this you are to determine
Q1) Does a shape1 parameter smaller than 1 
    make left tail go up or down.
Q2) Does a shape1 parameter large than 1 make the 
    make the left tail go up or down.

Q3) Does a shape2 parameter smaller than 1 
    make right tail go up or down.
Q4) Does a shape2 parameter larger than 1 make the 
    make the right tail go up or down.